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DETAILED ACTION 
Drawings 

1 . The drawings are objected to as failing to comply with 37 CFR 1 .84(p)(5) because they 
do not include the following reference sign(s) mentioned in the description on line 6 of page 7: 
33(2). Corrected drawing sheets in compliance with 37 CFR 1.121(d) are required in reply to the 
Office action to avoid abandonment of the application. Any amended replacement drawing sheet 
should include all of the figures appearing on the immediate prior version of the sheet, even if 
only one figure is being amended. Each drawing sheet submitted after the filing date of an 
application must be labeled in the top margin as either "Replacement Sheet" or "New Sheet" 
pursuant to 37 CFR 1 .121(d). If the changes are not accepted by the examiner, the applicant will 
be notified and informed of any required corrective action in the next Office action. The 
objection to the drawings will not be held in abeyance. 

Specification 

2. The lengthy specification has not been checked to the extent necessary to determine the 
presence of all possible minor errors. Applicant's cooperation is requested in correcting any 
errors of which applicant may become aware in the specification. 

Claim Objections 

3. Claims 1, 22, and 43 are objected to because of the following informalities: Claims 1, 22, 
and 43 contain the grammatical error "one a subdivision-inverse filter methodology or a least- 
squares optimization methodology is to be used." The office suggests the language "one of. ..." 
Appropriate correction is required. 
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Claim Rejections - 35 USC§ 103 

4. The following is a quotation of 35 U.S.C. 103(a) which forms the basis for all 

obviousness rejections set forth in this Office action: 

(a) A patent may not be obtained though the invention is not identically disclosed or described as set forth in 
section 102 of this title, if the differences between the subject matter sought to be patented and the prior art are 
such that the subject matter as a whole would have been obvious at the time the invention was made to a person 
having ordinary skill in the art to which said subject matter pertains. Patentability shall not be negatived by the 
manner in which the invention was made. 

5. Claims 1, 2, 9, 10, 16, 22, 23, 30, 31, 37, 43, 44, 51, 52, and 58 are rejected under 35 
U.S.C. 103(a) as being unpatentable over U.S. Patent No. 6,266,062 to Rivera (herein 
referred to as "Rivera") in view of Michael Loundsbery, Tony D. DeRose, and Joe Warren 
"Multiresolution Analysis for Surfaces of Arbitrary Topological Type," January 1997, 
ACM Transactions on Graphics, Vol. 16, No. 1, p. 34-73 (herein referred to as 
"Loundsbery et al") and in further view of Applicant's Admitted Prior Art: Denis Zorin, 
Peter Schroder, Wim Sweldens, "Interactive Multiresolution Mesh Editing," August 1997, 
Proceedings of the 24th Annual Conference on Computer Graphics and Interactive 
Techniques, p. 259-268 (herein referred to as "Zorin et al"). 

6. With regard to claim 22, Rivera et al discloses "a method of generating a coarse level 
mesh representation representing a surface, from a finer level mesh representation (abstract: 
"The dereflnement method, for each target vertex finding an associated set of neighbor vertices 
to be derefined; then eliminating each said vertex according an appropriate order such that the 
derefinement of said vertex allows to re-obtain a previous terminal edge whose bisection 
produced said vertex. ") including the steps of: 

a. indicator value generator step of, for respective ones of the points in the finer 
level mesh representation, evaluating an indicator function to generate an indicator value 
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(lines 18-32 of column 13: " The preferred embodiment of the derefinement method of 
this invention uses an integer-valued indicator function VE-IND (associated with the 
relation of precedence in the generation process between each vertex and its neighbors), 
such that, for each vertex VE, either (1) VE-IND is equal to zero whether VE is a vertex 
of the initial or improved mesh (not obtained by using the longest-edge refinement 
method of this invention), initialized in box 910 of FIG. 18... ") to determine a position for 
a corresponding point in the coarse level mesh representation; and 
b. a coarse level mesh generator step of determining, for each of the points that are 
to be provided in the coarse level mesh representation, a position in the coarse level mesh 
representation in response to the position of the corresponding point in the finer level 
mesh representation (lines 24-28 of column 16: 'Turning to FIG. 21, the derefinement 
method of this invention essentially comprises, for each target vertex VXto be derefmed 
or eliminated (box 1200) of associated node (VX, NX, GX) t where NX is the VE-IND 
value of VXand GX its generator edge (box 1220), ")...")> in accordance with the 
outcome of the indicator function (lines 37-40 of column 13: "The vertices whose VE- 
IND value is equal to 0 (and whose generator edge is equal to NULL) correspond to the 
initial mesh and cannot be derefined throughout the process. "). 
7. The indicator function in Rivera is analogous to the indicator function of claim 1 in that it 
is used to identify vertices resulting from a subdivision method (lines 23-25 of column 13: 
"...whether VE is a vertex of the initial or improved mesh (not obtained by using the longest- 
edge refinement method of this invention)... "). However, Rivera restricts the subdivision inverse 
step to meshes obtained from "longest-edge refinement," and does not teach a "subdivision 
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inverse filter." Loundsbery et al discloses a "subdivision inverse filter" (2 nd paragraph of section 
6.3: "Fowny multiresolution analysis the synthesis filters are defined by the relation ...and the 
analysis filters are obtained from the inverse relation...") to determine a position in the coarse 
level mesh representation in response to the position of the corresponding point in the finer level 
mesh representation (3 rd paragraph of section 6.3: "The analysis filters can be used to 
decompose a surface $ +1 (x) in given by ... into a lower resolution part in V/M 0 ) plus a 

detail part in W/M 0 ). "). 

8. Rivera and Loundsbery et al are analogous art because they are from a similar problem 
solving area: fine-to-coarse mesh generating methods. At the time of the invention, it would have 
been obvious to a person of ordinary skill in the art to use the method disclosed by Loundsbery et 
al to generate a coarse level mesh in the system disclosed by Rivera. The motivation for doing so 
would have been to reverse the subdivision methods other than the "longest edge refinement," 
such as those disclosed by Loundsbery et al in Table I on page 57. Therefore, it would have been 
obvious to combine Rivera with Loundsbery et al to obtain the invention specified in claim 22. 

9. However, neither Rivera nor Loundsbery et al disclose a "least squares optimization 
methodology." Rivera is unable to refine vertices that do not directly result from a subdivision 
method, as indicated by the indicator function (lines 37-40 of column 13: "The vertices whose 
VE-IND value is equal to 0 (and whose generator edge is equal to NULL) correspond to the 
initial mesh and cannot be derefined throughout the process. "). Zorin et al discloses using a least 
squares methodology in a fine-to-coarse level mesh generating method (2 nd paragraph of section 
3: "In this section we describe analysis which goes from fine to coarse. We first need smoothing, 
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i.e., a linear operation H to build a smooth coarse mesh at level /-/ from a fine mesh at level 
i... Least squares: One could define analysis to be optimal in the least squares sense... "). 

10. Rivera, Loundsbery et al and Zorin et al are analogous art because they are from a similar 
problem solving area: fine-to-coarse mesh generating methods. At the time of the invention, it 
would have been obvious to a person of ordinary skill in the art to incorporate a least squares 
optimal method as disclosed by Zorin et al to refine the vertices that are indicated by the 
indicator function to correspond to the initial mesh {lines 37-40 of column 13) in the system 
disclosed by Rivera. The motivation for doing so would have been to generate a coarser mesh 
that would otherwise possible under the condition that only vertices that are the direct result of 
the subdivision can be derefined. Therefore, it would have been obvious to combine Zorin et al 
with Rivera and Loundsbery et al to obtain the invention specified in claim 22. 

1 1 . With regard to claim 23, Zorin et al further discloses "Laplacian generator step of 
generating a Laplacian value for said respective ones of the points in the finer level mesh 
representation" (2 nd paragraph of section 3: "Use the average... to define the discrete Laplacian 
H:= (I + ju£)(I+A,L); 4 th paragraph of section 3: "As indicated in Fig. 7 the detail vectors are 
defined as (J -Ss*' 1 *)=(F*) { f (I-SH)sf, i.e. the detail vectors at level i record how much the 
points at level i differ from the result of subdividing the points at level /-/. ")■ 

12. With regard to claims 30, 3 1 , and 37, Zorin et al further discloses "the coarse level mesh 
generator step includes the step of determining, for at least one of the points that are to be 
provided in the coarse level mesh representation, the position in the coarse level mesh 
representation as the position of the corresponding point in the finer level mesh representation 
(6 th paragraph of section 4: "Analysis computes points on the coarser level i-1 using smoothing 
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(smooth), subdivides /' (subd), and computes the detail vectors d (cf. Fig. 7). ") if the 
magnitude of the Laplacian value generated during the Laplacian generator step is below a 
predetermined threshold value" (1 st paragraph of section 4.1: "Adaptive analysis avoids the 
storage cost associated with detail vectors below some threshold e A by observing that small 
detail vectors imply that the finer level almost coincides with the subdivided coarser level. The 
storage savings are realized through tree pruning. 2 nd paragraph of section 4. 1 : 'Triangles 
that do not contain details above the threshold are unrefined,.. "). 

13. With regard to claims 30, 31, and 37, Zorin et al does not chose between a subdivision 
inverse filter or least squares methodology when the detail vectors are below a given threshold. 
However, one of ordinary skill in the art would recognize the magnitude of the detail vector 
indicates which method would be more suitable from the statement in the 1 st paragraph of section 
4.1 : "small detail vectors imply that the finer level almost coincides with the subdivided coarser 
level." At the time of the invention, it would have been obvious to a person of ordinary skill in 
the art to further modify the combination of Rivera, Loundsbery et al, and Zorin et al to use a 
Laplacian as disclosed by Zorin et al to generate a detail vector to be used to determine the type 
of coarse level generation step, and use the subdivision inverse filter in the instance the detail 
vector indicates the fine mesh is subdivided coarser level, and the least squares in the instance it 
does not. The motivation for doing so would have been to enable the indicator function of Rivera 
to determine finer levels coinciding with coarse levels for subdivision schemes other than those 
generated by "longest-edge refinement " Therefore, it would have been obvious to further 
modify the combination of Rivera, Loundsbery et al, and Zorin et al to obtain the invention 
specified in claims 23, 30, 31, and 37. 
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14. Claim 1 is rejected with the rationale of claim 22. Claim 1 is similar in scope to claim 22. 
Rivera shows a fine-to-coarse level mesh generating arrangement (line 6 to column 14 through 
lines 4 of column 4: " The actualization of the longest-edge mesh data structure, either by point 
insertion (refinement) or by point elimination (derefmement) is locally performed by the parallel 
computers of FIG. 20 and FIG. 23 respectively which in turn modify the shared mesh data 
structure according the diagrams of FIGS. 19 and 24, respectively. "). 

15. Claims 2, 9, 10, and 16 are rejected with the rationale of claims 23, 30, 31, and 37, 
respectively. Claims 2, 9, 10, and 16 are similar in scope to claims 23, 30, 31, and 37. 

1 6. Claim 43 is rejected with the rationale of claim 22. Claim 43 is similar in scope to claim 
22. The feature of a computer program product is inherent in the system disclosed by Rivera in 
Figure 20, as the system would be inoperable without a computer program product instructing 
the computer to manipulate the data structure in the manner defined by the method disclosed by 
Rivera. 

1 7. Claims 44, 5 1 , 52, and 58 are rejected with the rationale of claims 23, 30, 3 1 , and 37, 
respectively. Claims 44, 51, 52, and 58 are similar in scope to claims 23, 30, 31, and 37. 

Allowable Subject Matter 

18. Claims 3-8, 11-15, 17-21, 24-29, 32-36, 38-42, 45-50, 53-57, and 59-63 are objected to as 
being dependent upon a rejected base claim, but would be allowable if rewritten in independent 
form including all of the limitations of the base claim and any intervening claims. 

1 9. The following is a statement of reasons for the indication of allowable subject matter: 
While the prior art of Rivera, Loundsbery et al, and Zorin et al show the limitations recited in the 
parent claims, the mathematical expressions claimed in 3-7, 11-15, 17-21, 24-28, 32-36, 38-42, 
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45-49, 53-57, and 59-63 are not shown nor are they obvious in view of the teachings of the prior 
art. 

Conclusion 

20. The prior art made of record and not relied upon is considered pertinent to applicant's 
disclosure. U.S. Patent No. 5,506,947 to Taubin et al discloses using a Gaussian kernel in a 
surface smoothing method. U.S. Patent No. 6,009,435 to Taubin et al discloses a multiresolution 
polygonal mesh. 

Any inquiry concerning this communication or earlier communications from the 
examiner should be directed to Jason M. Repko whose telephone number is 571-272-8624. The 
examiner can normally be reached on Monday through Friday 8:30 am -5:00 pm. 

If attempts to reach the examiner by telephone are unsuccessful, the examiner's 
supervisor, Ulka Chauhan can be reached on 571-272-7782. The fax phone number for the 
organization where this application or proceeding is assigned is 571-273-8300. 

Information regarding the status of an application may be obtained from the Patent 
Application Information Retrieval (PAIR) system. Status information for published applications 
may be obtained from either Private PAIR or Public PAIR. Status information for unpublished 
applications is available through Private PAIR only. For more information about the PAIR 
system, see http://pair-direct.uspto.gov. Should you have questions on access to the Private PAIR 
system, contact the Electronic Business Center (EBC) at 866-217-9197 (toll-free). 



JMR 



ULKA CHAUHAN 



SUPERVISORY PATENT EXAMINER 
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Multiresolution Analysis for Surfaces of 
Arbitrary Topological Type 

MICHAEL LOUNSBERY, TONY D. DeROSE, and JOE WARREN 
University of Washington, Seattle 



Multiresolution analysis and wavelets provide useful and efficient tools for representing 
functions at multiple levels of detail. Wavelet representations have been used in a broad range 
of applications, including image compression, physical simulation, and numerical analysis. In 
this article, we present a new class of wavelets, based on subdivision surfaces, that radically 
extends the class of representable functions. Whereas previous two-dimensional methods were 
restricted to functions denned on R 2 , the subdivision wavelets developed here may be applied 
to functions denned on compact surfaces of arbitrary topological type. We envision many 
applications of this work, including continuous level-of-detail control for graphics rendering, 
compression of geometric models, and acceleration of global illumination algorithms. Level-of- 
detail control for spherical domains is illustrated using two examples: shape approximation of 
a polyhedral model, and color approximation of global terrain data. 

Categories and Subject Descriptors: G.1.2 [Approximation]: Spline Approximation; 1.3.5 
[Computer Graphics]: Computational Geometry and Object Modeling — surfaces and object 
representations; J. 6 [Computer-Aided Engineering]: Computer-Aided Design (CAD) 

General Terms: Algorithms, Design, Theory 

Additional Key Words and Phrases: Compression, geometric modeling, level-of-detail control, 
splines, subdivision surfaces, wavelets. 



* 1. INTRODUCTION 

Multiresolution analysis and wavelets have received considerable attention 
in recent years, fueled largely by the diverse collection of problems that 
benefit from their use. The basic idea behind multiresolution analysis is to 
decompose a complicated function into a "simpler" low resolution part, 
together with a collection of perturbations, called wavelet coefficients, 
necessary to recover the original function. For many of the functions 
encountered in practice, a large percentage of the wavelet coefficients are 
small, meaning that good approximations can be obtained by using only a 
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few of the largest coefficients. Impressive "lossy* compression rates for 
images have been achieved using this type of approximation [DeVore et al. 

1992]. . 

There are many constructions of wavelets for functions parametrized 
over the interval [Andersson et al. 1993; Chui and Quak 1992; Cohen et al. 
1993; Meyer 1992]. These have found use in signal processing [Mallat 
1989], B-spline modeling [Finkelstein and Salesin 1994], motion planning 
[Liu et al. 1994], and many other applications involving functions parame- 
trized in one dimension. 

Two-dimensional wavelets are important for a variety of applications 
including image compression. They are generally constructed by forming 
tensor products of univariate wavelets [Daubechies 1992], in much the 
same way that tensor product B-spline surfaces are formed by products of 
univariate B-splines. Unfortunately, tensor-product constructions require 
that the functions to be decomposed be defined on U 2 or on a periodic 
version of R 2 , that is, the cylinder or torus. There also exist nontensor- 
product constructions for wavelets on R 2 [Daubechies 1992; Jia and Micch- 
elli 1991], but none of these methods are applicable to functions defined on 
more general topological domains, such as spheres. Thus, existing methods 
are not well suited for decomposing and compressing surfaces such as the 
ones shown in Figures 10 and 12, since they are described by parametric 
functions on the sphere. 

In this article we show that by using techniques from subdivision 
surfaces, multiresolution analysis can be extended to functions defined on 
domains of arbitrary topological type (the topological type of a two-dimen- 
sional surface or domain refers to its genus, presence of boundary curves, 
etc.). This generalization, which we term "subdivision wavelets," dramati- 
cally extends the class of applications to which multiresolution analysis can 
be applied, including: 

—Polyhedral compression. Using wavelet-based techniques, most polyhe- 
dral models may be compressed to yield a more compact approximation. 
Compression saves both storage space and the time that is required to 
process a surface model. This article develops efficient wavelet compres- 
sion methods for surfaces similar to those that have proven effective for 
images and one-dimensional functions. These techniques are capable of 
L°° and L 2 approximation, and run more quickly than many other 
compression methods. This application is explored in more detail in Eck 
et al. [1995]. 

—Continuous level-of -detail control. When a complex shape is rendered in 
an animation, a fully detailed representation of the shape contains much 
more detail than is required for all but the closest view. Using a 
compressed subdivision wavelet representation of complex objects, it is 
possible to greatly reduce the number of polygons in a scene without 
significantly impacting the visible detail (see Figure 10). Moreover, it is 
possible to smoothly vary the detail, avoiding the discontinuous jumps 
that occur when suddenly switching between distinct models. This appli- 
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cation is discussed in more detail in Section 7.5. Nonwavelet treatments 
of this problem may be found in Turk [1992], Schroeder et al. [1992], and 
Hoppe et al. [1993]. 

— Compression of functions defined on surfaces. Consider the situation 
shown in Figure 12 where a globe (a geometric sphere) is pseudo-colored 
according to elevation. The pseudo-coloring can be thought of as a 
function that maps each point on the sphere to an RGB triple. A 
straightforward method for storing the function is to store its value at a 
large number of regularly distributed points; in this case more than one 
million points were used. The methods in this article can be used to 
create compressed wavelet approximations of varying complexity. (The 
mesh lines on the original surface are so dense that the image shown in 
Figure 11(b) is nearly black.) 

—Multiresolution editing of surfaces. Hierarchical B-splines, as introduced 
by Forsey and Bartels [1988], provide a powerful mechanism for editing 
shapes at various levels of detail. However, hierarchical B-splines can 
only represent a restricted class of surface topologies. The methods 
described here provide an alternative to hierarchical B-splines, and are 
capable of representing smooth multiresolution surfaces of arbitrary 
topological type. Editing at fractional levels of detail can also be achieved 
using the methods developed by Finkelstein and Salesin [1994]. 

—Surface optimization. The multiple levels of approximation produced by 
wavelet techniques offer a sort of multigrid technique for optimization. 
Pentland [1992] uses wavelet methods to implement multigrid optimiza- 
tion for surface interpolation over regular grids. Meyers [1994a; 1994b] 
shows how wavelets can accelerate the reconstruction of surfaces from 
contour data. This previous work suggests that subdivision wavelets may 
find use in optimization over surfaces of arbitrary topological type. 

— Numerical solution of integral and differential equations. Wavelet repre- 
sentations of functions on surfaces of arbitrary topological type appear 
well suited for adaptive numerical solution, along the lines of Beylkin et 
al. [1991]. 

A carefully constructed special-purpose algorithm may produce superior 
results for some of these applications. However, the above examples indi- 
cate the versatility of subdivision wavelets, and show them to be a reusable 
tool that naturally and efficiently solve or accelerate a wide range of 
common problems in computer graphics and modeling. As an illustration of 
this versatility, a complex object of arbitrary topological type may be 
modeled at multiple levels of detail, stored in a compressed form, viewed in 
an animation, and rendered using global illumination — all naturally imple- 
mented using the same underlying surface wavelet representation. 

This article presents a theoretical foundation for developing multiresolu- 
tion analysis for surfaces of arbitrary topological type. (Additional details of 
implementation and an expanded treatment of applications may be found 
in Lounsbery [1994], Eck et al. [1995], and Certain et al. [1996].) 
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Wavelet coefficients Wavelet coefficients 

Fig. 1. Decomposition of a polyhedral surface. 



The construction of subdivision wavelets described herein applies directly 
to a variety of existing subdivision schemes. These include piecewise linear 
subdivision (producing polyhedra); the schemes of Catmull and Clark 
[1978], Halstead et al. [1993], Loop [1987] or Dyn et al. [1990] (producing 
tangent-plane-smooth G 1 subdivision surfaces); and the modifications by 
Hoppe [1994] and Hoppe et al. [1994] (producing piecewise smooth surfaces 
with selected discontinuities). More generally, the techniques presented 
here may be used to construct wavelets for any local, stationary, continu- 
ous, uniformly convergent subdivision scheme (these terms will be de- 
scribed in Section 4). 

The remainder of the article is structured as follows. In Section 2, we 
present a high-level preview of how to convert a polyhedral object to 
multiresolution form. In Section 3, we provide some background on mul- 
tiresolution analysis. In Section 4, we provide a brief overview of subdivi- 
sion surfaces, and we show that they can be used to define a collection of 
scaling functions necessary for multiresolution analysis. In Section 5, we 
show that inner products of these scaling functions can be computed 
exactly. In Section 6, we use the inner products to construct wavelets, and 
to construct locally supported approximations to them. In Section 7, we 
apply the theory to polyhedral compression. In Section 8, we apply the 
theory to the problem of compression and editing of the smooth surfaces 
described by Dyn et al. [1990], Finally, in Section 9, we summarize the 
contributions of the article, and list several topics for future work. 

2. A PREVIEW OF THE METHOD 

Although the mathematical underpinnings of multiresolution analysis of 
surfaces are somewhat involved, the resulting algorithms are quite simple. 
Before diving into the details, we give here a brief description of how the 
method can be applied to decompose the polyhedral object shown in Figure 
Ka). 

As mentioned in Section 1, a main idea behind multiresolution analysis is 
the decomposition of a function (in this case a polyhedron expressed as a 
parametric function on the sphere) into a low resolution part and a "detail" 
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part. The low resolution part of the polyhedron in Figure 1(a) is shown in 
Figure Kb). The vertices in Figure 1(b) are computed as certain weighted 
averages of the vertices in Figure 1(a). These weighted averages essentially 
implement a low pass filter denoted as A. The detail part consists of a 
collection of fairly abstract coefficients, called wavelet coefficients, com- 
puted as weighted differences of the vertices in Figure 1(a). These differ- 
encing weights form a high-pass filter B. The decomposition process, 
technically called analysis, can be used to further split Figure 1(b) into an 
even lower resolution version and corresponding wavelet coefficients. This 
cascade of analysis steps, referred to as a filter bank algorithm, culminates 
with the coarsest-level representation in Figure 1(c), together with wavelet 
coefficients at each level. 

The analysis filters A and B are constructed so that the original polyhe- 
dron can be recovered exactly from the low-resolution version and the 
wavelet coefficients. Recovery, technically called synthesis, reconstructs 
Figure 1(a) from Figure Kb) and the finest-level wavelet coefficients. 
Recovery refines each triangle of Figure Kb) into four subtriangles by 
introducing new vertices at edge midpoints, followed by perturbing the 
resulting collection of vertices according to the wavelet coefficients. The 
refining and perturbing steps are described by two other filters P (the 
refining filter) and Q (the perturbing filter), collectively called synthesis 
filters. 

The trick is to develop the four analysis and synthesis filters so that: (1) 
the low-resolution versions are good approximations of the original object 
(in a least-squares sense); (2) the magnitude of a wavelet coefficient reflects 
a coefficient's importance by measuring the error introduced when the 
coefficient is set to zero; and (3) analysis and synthesis filter banks should 
have time complexity that is linear in the number of vertices. 

3. BACKGROUND ON MULTIRESOLUTION ANALYSIS 

Multiresolution analysis as formulated by Mallat [1989] and Meyer [1993] 
provides a convenient framework for developing the analysis and synthesis 
filters. There are two basic ingredients for a multiresolution analysis: an 
infinite chain of nested linear function spaces V° C V 1 C V 2 C * * • and an 
inner product if, g) defined on any pair of functions /,gG V 7 , for some j < 
oo. Intuitively, V J contains functions of resolution./, with the detail increas- 
ing as j increases. 

The inner product is used to define the orthogonal complement spaces W J 
as 

W := {/E V* +l \<f 9 g) = 0 g E V}. 

Orthogonal complements are often written as V J+1 = V J © W J because any 
function f J+1 E V J+l can be written uniquely as an orthogonal decomposi- 
tion 
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where f J E V J and h J E W J . Orthogonal decompositions are important for 
approximation purposes: it is easy to show that/ 7 is the best approximation 
to / v+1 in that it is the unique function in V J that minimizes the least- 
squares residual (f J+1 - f J f f J+1 - f J ). Thus, given a high-resolution 
function f J+1 , its low-resolution part is f J , and its detail part is W . 

The following terminology is now standard: scaling functions refer to 
bases for the spaces V j y and wavelets refer to bases for the orthogonal 
complement spaces. As shown in Section 6.3, the analysis and synthesis 
filters are determined by considering various ways of changing bases 
between scaling functions and wavelets. 

The strictest form of wavelets, such as those constructed by Daubechies 
[1988] and Mallat [1989], are fully orthogonal, meaning that every wavelet 
i^(x) is orthogonal to every other wavelet i/r£(x); that is, <iK( x )» iM'( x )> = 0 
unless j = j' and i = i'. 

It has been shown to be impossible to construct wavelets that are 
simultaneously locally supported, fully orthogonal, and symmetric. It is 
therefore sometimes convenient to relax the fully orthogonal condition to 
the semiorthogonal setting, and require only orthogonality between wave- 
lets at different levels: <*M(x), i/rf,(x)) = 0 unless j = f. The locally 
supported B-spline wavelets constructed by Chui [1992] are a good example 
of a semiorthogonal construction. 

Finally, the least restrictive form of wavelets are said to be biorthogonal, 
a term coined by Cohen et al. [1992] to refer to the setting where 
orthogonality is dropped entirely. It is not even necessary for the the 
wavelet spaces W to be orthogonal complements — in general W J is just 
some complement of V 7 in V J+1 . Our construction, described in Section 6, 
results in locally supported biorthogonal wavelets. 

4. NESTED LINEAR SPACES THROUGH SUBDIVISION 

A fundamental requirement for multiresolution analysis is a sequence of 
nested linear spaces. In this section, we carry this property to surfaces of 
arbitrary topological type, demonstrating the existence of scaling functions 
on subdivision surfaces. 

The nested sequence of linear spaces required by multiresolution analy- 
sis are ordinarily obtained by defining a single scaling function <f>(x) that 
satisfies a refinement equation of the form 

*(*) = 2 Pi*(2x - i) 

for some fixed constants p t . The refinement equation (sometimes called a 
two-scale relation) guarantees that the spaces defined as 

V s : = Span{<M2'x - = -«>, - - . , 

are nested. In other words, the nested spaces are generated by translations 
and dilations of a single refinable function <£>(jc). 
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To generalize these ideas to domains of arbitrary topological type, one 
could attempt to make definitions for what it means to dilate and translate 
a function on an arbitrary topological domain. One could then try to find a 
refinable scaling function and proceed as before to define orthogonal 
complements, wavelets, and so on. We have instead chosen what appears to 
be a simpler approach. 

In this section, recursive subdivision is used to define a collection of 
functions <f>i(x) that are refinable — in the sense that each function with 
superscript j lies in the span of the functions with superscript j + 1; the 
argument x is a point that ranges over a domain 2-manifold of arbitrary 
topological type. In one respect, this is a generalization of the approach 
taken by Daubechies [1992], whose locally supported orthogonal scaling 
functions are also defined through a recursive subdivision procedure. 
Although in general the <f> J+1 (x) are not simple dilates of the ^(x), one can 
nonetheless use them to define a sequence of nested spaces. 

(We should note that the realization that translation and dilation are not 
strictly necessary for the construction of wavelets has also been noted 
independently by Sweldens [1994] and Carnicer et al. [1994].) 

4.1 Subdivision Surfaces 

Intuitively speaking, subdivision surfaces are defined by iteratively refin- 
ing a control polyhedron M° so that the sequence of increasingly faceted 
polyhedra M 1 , M 2 , . . . converge to some limit surface M 00 . In each 
subdivision step, the vertices of M* 7 * 1 are computed as affine combinations 
of the vertices of M j . Thus, if V is a matrix whose ith row consists of the x, 
y, and z coordinates of vertex i of M J , there exists a nonsquare matrix of 
constants P- 7 such that 

V'" +1 = F'V. (1) 

The matrix P J therefore characterizes the subdivision method. The beauty 
of subdivision surface schemes is that the entries of P J depend only on the 
connectivity of the vertices in M\ not on the geometric positions of the 
vertices. Subdivision schemes are typically local, meaning that each vertex 
of M J+1 is computed as an affine combination of nearby vertices of M J . 

The simplest example of such a scheme is polyhedral subdivision. Given a 
polyhedron M° with triangular faces, a new polyhedron M 1 is built by 
splitting each triangular face of M° into four subfaces as in Figure 2. The 
matrix P° characterizing the first subdivision step is also shown in Figure 
2. Running this subdivision scheme for j steps on an initial triangular mesh 
M° produces a mesh M J . M J includes the vertices of M 0 together with new 
vertices introduced through subdivision. The valence (the number of edges 
incident to a vertex) of the vertices of M J corresponding to the original 
vertices in M° remains fixed. The new vertices introduced through subdivi- 
sion however are always of valence six, corresponding to a regular triangu- 
lar tiling of the surface. As the mesh is further subdivided, the so-called 
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Fig. 2. Polyhedral subdivision of a tetrahedron. The matrices shown here represent various 
filters associated with the scheme, and are explained in Section 6. 



extraordinary points (any original vertex of valence other than six) become 
increasingly isolated in an otherwise regular tiling of the surface. 

Polyhedral subdivision converges to the original polyhedron covering M°, 
that is, to a C° surface. However, other schemes have been developed that 
converge to tangent-plane smooth limit surfaces that either approximate or 
interpolate the vertices of M°. Subdivision schemes can be further catego- 
rized as being either primal or dual. A subdivision scheme is primal if the 
faces of the mesh are split into subfaces by the refinement procedure. 
Catmull and Clark subdivision [1978; Halstead et al. 1993] is a primal 
scheme based on subdivision of quadrilateral faces. Polyhedral subdivision, 
the w butterfly ,, scheme of Dyn et al. [1990] and Loop's method [1987] are 
primal schemes based on subdivision of triangular faces. A scheme is dual 
if the structure of the refined mesh is obtained by doing a primal subdivi- 
sion followed by taking the polyhedral dual of the result. Doo and Sabin 
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Fig. 3. The subdivision limiting process. 



[1978] and Doo [1978] subdivision is a dual scheme based on quadrilaterals. 
Additionally, a stationary subdivision scheme is one in which there exists 
an integer j such that the same subdivision rule is used for refinement 
steps numbered larger than j. 

Although one can certainly construct fractal-like subdivision schemes 
that do not converge to well defined surfaces, in this article we shall 
consider only local schemes that uniformly converge to a continuous surface 
no matter where the control vertices are placed. In the remainder of this 
article, subdivision schemes will be assumed to be local, stationary, contin- 
uous, and uniformly convergent. For simplicity we shall restrict the discus- 
sion to primal triangular subdivision schemes, although our results also 
hold for primal quadrilateral schemes. 

4.2 Refinable Scaling Functions Through Subdivision 

We show here that subdivision can be used to define a collection of 
refinable scaling functions. We do this by first showing that subdivision 
surfaces can be parametrized by a function S(x), where x ranges over M°. 
We then show that the parametrization can be used to define the scaling 
functions. Parametrizing the scaling functions over a domain M 0 of arbi- 
trarily topological type differs sharply from the more usual method of 
parametrizing surface basis functions over a piece of the plane. 

In general terms, a surface parametrization is nothing more than a 
correspondence between points in a two-dimensional domain and points on 
the surface. The idea behind establishing a parametrization for a subdivi- 
sion surface is to track a point x on M 0 through the subdivision process (see 
Figure 3). The point x being tracked will converge to a point on the limit 
surface, thereby establishing a correspondence S between points on M° and 
points on the limit surface. We now make these ideas more precise. 

For primal subdivision schemes, the refinement step that carries the 
mesh M 5 " 1 into M 8 consists of two substeps: the splitting step and the 
averaging step. In the splitting step, each face of M 5 " 1 is split into four 
subfaces by introducing vertices at midpoints of edges, creating auxiliary 
mesh M 8 , as shown in Figure 3. The averaging step uses local weighted 
averaging to compute the vertex positions of M 8 from the vertex positions of 
M 8 . All primal subdivision schemes share the splitting step — they differ 
only in the weights used in the averaging step. 
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The parametrization S(x) can now be established by using a limiting 
process. The limiting process producing Six) consists of three steps: 

(1) S°(x) := x,x GM°. 

(2) Suppose thatS fi_1 (x) lies in triangle (t3*, v% 9 v%) oiSl 8 with barycentric 
coordinates (a, )3, 7). Then 

S*(x) :=a^-h^J + 7^, (2) 
where (v 8 ay v%, v 8 c ) is the triangle of M 8 corresponding to (0 8 a , 6 8 b , 0 S C ) of 

iff*. 

(3) S(x) := hm^S^xl 

Next, we show that the parametrization S(x) induces a collection of 
refinable scaling functions. 

Lemma 4.2.1. For all j ^ 0, all s ^ j, and all i ranging over the vertices 
ofM J , there exist functions <fi- J : M° -» R such that 

s*(x) = s^i+rw. 

i 

Proof. It is convenient to write the statement of the lemma in matrix 
form as 

S s (x) = **^"(x)V\ 

where $> s ~~ J (x) is the row vector whose ith component is <#(x). 

The linear combination of Equation 2 above can be rewritten in matrix 
notation as 

S*(x) = b a (x)V s , 

where b*(x) is the barycentric coordinate vector of x with respect to M 8 \ 
that is, 

b*(x) = (0 • • • OaO • • • 000 • • • OyO ♦ • * 0), 

where a occurs at index a, j3 occurs at index 6, y occurs at index c. At each 
refinement step k = 1, . . . , s, the vertices of M k can be computed from 
affine combinations of the vertices of M k ~ x . Therefore, there must exist a 
chain of (non-square) matrices P°, . . . , P* -1 such that 

V* = p*-ip*-2 . . . p°v°. 

Thus, 

S*(x) = b^xJP 5 - 1 ? 5 - 2 • • • P'V. 
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The desired result follows immediately by making the definition 

d> a ^'(x) : =b 5 (x)P al P 5 - 2 • • • F>. □ 

As a simple corollary to Lemma 1, we note that 

4>— >( x ) = 4>^ +1 (x)P'. (3) 

Theorem 4.2.1. For any subdivision procedure, and for any j s 0, there 
eocist scalar-valued functions <t>{(x), x e such that 

i 

Proof. Using Lemma 1 and the definition of Six): 

S(x) = bm(**^V0. 

Since the subdivision scheme is assumed to be local and uniformly conver- 
gent for any choice of control points, S(x) must exist if we choose all entries 
of V to be the origin, except for v{ 9 the ith entry of V. This choice of control 
points leads to a surface of the form 



lim <f>r'(xM = v*A lim <f>r y (*) 

Since the surface is well defined, the sequence on the right must converge 
to u-^Cx), where 

<M(x) := lim <J>r'(x). 

By linearity, the surface for an arbitrary set of control points can be written 
as in Equation 4. □ 

We will find it convenient to rewrite Equation 4 in matrix form as 

S(x) = <Hx)V, (5) 

where & J (±) denotes the row matrix of scaling functions <M(x), and where 
V is as in Equation 1. Equation 5 shows an analogy with B-splines, where 
the <H(x) are comparable to the B-spline basis functions, and the V are 
akin to the control points. 

As a corollary to Theorem 4.2.1, if the subdivision scheme generates 
continuous surfaces, the scaling functions <f>{ are continuous, hence they 
are also integrable [Bartle 1964]. 

We may now establish the refinability of the scaling functions defined in 
Theorem 4.2.1. 

ACM Transactions on Graphics, Vol. 16, No. 1, January 1997. 



Multiresolution Analysis • 45 



Theorem 4.2.2. The scaling functions <H(x) are refinable. 

Proof. Starting with Equation 3 and taking limits as s tends toward 
infinity, it follows from the existence of the <H(x) that 

*'"(x) = 4>'' +1 (x)PA (6) 

This equation establishes refinability since it states that each of the 
functions <t>{(x) can be written as a linear combination of the functions 

<M +1 (x). □ 

For primal subdivision schemes, it is convenient to write Equation 6 in 
block matrix form by writing 4>- ;+1 (x) as 

*> +1 (x) = (0 ;+1 (x) >T> +1 (x)), (7) 

where 0- /+1 (x) consists of all scaling functions <H +1 (x) associated with the 
old vertices of M J (the black vertices in Figure 2) and J( J+1 (x) consists of 
the remaining scaling functions associated with the new vertices added 
when obtaining M J+1 from M j (the white vertices in Figure 2). Equation 6 
can now be expressed in block matrix form: 

*>(x) = Jf> +1 (x))(j},), (8) 

where O- 7 and N J represent the portions of the subdivision matrix P 7 which 
weight the "old" and "new" vertices, respectively. The block matrix decom- 
position of P° for the example tetrahedron appears in Figure 2. 

4.3 Nested Linear Spaces 

Given these relations, a chain of nested linear spaces V J (M°) associated 
with a mesh M 0 can now be defined as follows: 

V(M°) : = Span(4>>(x)), 

where the V J (M°) are spaces of scalar-valued functions. 

Equation 6 implies that these spaces are indeed nested; that is, 

V°(M°) C V l (M°) C • - • . 

The notation V J (M°) is to emphasize that the linear spaces are adapted to 
M° in that they consist of functions having M° as the domain. 

5. INNER PRODUCTS 

Given a chain of nested linear spaces, the other necessary ingredient for 
the creation of a multiresolution analysis is the existence of an inner 
product on these spaces. In this section, we define an inner product and 
give a method for exactly computing the inner product of any pair of 
functions defined through subdivision. 
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Inner products between pairs of scaling functions are used in the con- 
struction of wavelets, as described in Section 6. For an efficient implemen- 
tation of subdivision wavelets, neither the inner products nor the wavelets 
need to be computed at run time, if the base mesh M° is known in advance. 

5.1 Definition 

Let two functions f, g E V J (M°) 9 j < 00 be linear combinations of 
(scalar-valued) scaling functions defined through subdivision, as in Section 
4. We define the inner product of / and g to be 



where the area form dx is taken to be the area for a triangulation 
homeomorphic to M° consisting of equilateral triangles with equal area. 
Equivalently, the inner product can be expressed as 



where A(M°) denotes the set of triangular faces of M°, and where dx' is the 
usual Euclidean area form for the triangle t in R 3 . 

Our definition of inner product implies that triangles of different geomet- 
ric size and shape are weighted equally; that is, the inner product is 
independent of the geometric positions of the vertices of M°. This inner 
product has two important consequences. 

First, in the process of constructing the least-squares best wavelet 
approximation to a function, each approximated triangle is weighted 
equally, independent of its geometric size. Although this simplification 
ignores the fact that for most data sets all triangles are not really the same 
size, we have obtained good results for many examples, including those 
described in Section 7.5. 

Second, because all triangles can be treated equally, the wavelet spaces 
are invariant of the geometry of the mesh and filters only change at the 
extraordinary points. Thus, a significant amount of precomputation of 
inner products and wavelets can be performed, allowing the wavelet 
algorithms to be implemented more efficiently. 

An alternative is to define the inner product so as to weight the integral 
by the areas of triangles in M°. Whether such a definition has enough 
important practical benefit to offset its much increased computation may be 
an interesting topic for future research. A step in this direction has already 
been made by Schroder and Sweldens [1995]. 

5.2 Computation 

For piecewise linear subdivision (leading to polyhedral surfaces), the scal- 
ing functions <f>j(x) are simply the hat functions over M°. For the case when 
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functions f and g are combinations of piecewise linear scaling functions, it 
is therefore fairly simple to directly compute the integral of Equation 9. 

In general, however, the lack of a closed form for the scaling functions 
makes explicit integration quite difficult. In these cases, one could estimate 
the inner product (f 9 g) by subdividing the scaling functions some number 
of times and directly integrating the resulting piecewise linear approxima- 
tion. In this section, we will see that such estimation is unnecessary, and 
that exact integration is still possible for the general case. 

One way to compute inner products is to follow the approach in Halstead 
et al. [1993]. Their approach for computing integrals over subdivision 
surfaces was to observe that away from extraordinary points, integrals 
could be computed exactly since Catmull-Clark subdivision converges to 
uniform bicubic B-splines in regular regions (those regions removed from 
extraordinary points). They also observed that the subdivision process 
could then be used to compute an integral over the entire surface by 
evaluating a geometric series. To generalize their method to other schemes, 
one must first compute inner products in regular regions. This is relatively 
easy for some schemes, such as Loop's, that converge to polynomials in 
regular regions, but it is harder for schemes such as the butterfly method 
that are nowhere polynomial. 

A general method for integration in regular regions was recently devel- 
oped by Dahmen and Micchelli [1993]. It operates by reducing the problem 
of computing inner products to one of computing eigenvectors of a matrix 
defined by the refinement equations. 

Although it is possible to combine the methods of Halstead et al. [1993] 
and Dahmen and Micchelli [1993], we present here a method that is 
somewhat simpler, in that it computes inner products directly as the 
solution to a homogeneous system of linear equations. 

Let f and g be functions given as expansions in <b\\ 



Bilinearity of the inner product allows (f, g) to be written in matrix form as 

if, 8) = gW, 

where f and g are column matrices consisting of the coefficients off and g> 
respectively, and where V is the square matrix whose i, t'-th entry is {\ J \ v — 

The ith row of V contains the inner product of with each of the other 
scaling functions 4>i'> It is convenient to view these entries geometrically by 
constructing an inner product mask around each vertex. The inner product 
mask for the ith vertex of M j assigns to each vertex V of M J a multiple of 
the scalar (<f> j i9 <f> J v ). 

In the case of polyhedral subdivision, explicit calculation leads to the 
inner product mask around a vertex of valence n shown in Figure 4, where 
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Fig. 4. The polyhedral inner product mask around a vertex of degree n. 



the central weight is simply the valence of the central vertex. As an 
example, the complete inner product matrix 1° for the tetrahedron appears 
in Figure 2. Note that in this case the mask entries must be divided by 3 to 
obtain the matrix entries. 

For more general subdivision procedures, such as the butterfly scheme, 
the limit surface has no closed form, precluding a brute-force explicit 
integration. However, as we now show, it is possible to exactly determine 
the entries of P by solving a linear system, without resorting to numerical 
integration. We will first present the general integration procedure, and 
then illustrate it using a simple one-dimensional example. 

The key to the linear system is to observe that when all triangles at 
subdivision level j are weighted equally, as was done in the inner product 
formulation in Section 5.1, a recurrence relation exists between P and P +1 . 
Specifically, P can be written as 



where the integrand represents a matrix outer product, and where the 
integral of a matrix of functions is defined to be the matrix of integrals. The 
refinement property of Equation 6 can now be used to establish the 
recurrence 



Since the subdivision rules are local, the corresponding scaling functions 
are locally supported, and the support of a scaling function <H(x) shrinks as 
j increases. Thus, beyond some level f the support of each of the scaling 
functions <f>{ will contain at most one extraordinary point. Furthermore, 
the local neighborhood at level j' + 1 will be a shrunken version of the 
local neighborhood at level j'. Hence, the scaling functions in 4> J '(x) and 
4> J +1 (x) are closely related. 
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Fig. 5. The diagram on the left depicts the triangulation of M° produced after J' subdivision 
steps, that is, after j' recursive midpoint splits; the diagram on the right corresponds toy' + 1 
subdivision steps. denotes the scaling function for the i-th vertex of M J \ and <#-'. +1 denotes 
the scaling function for the corresponding i'-th vertex of Af*' + 1 . The map a t is such that the 
bary centric coordinates of x and o^x) within their respective surrounding triangles are equal. 



More precisely, for each scaling function <f>{ there is at least one V such 
that 

<M'(x) = Vi(x)), 

where a t : M° -> M° is a piecewise linear function that maps triangles in 
the support of into the corresponding triangles of <f>{> + 1 9 as depicted in 
Figure 5. As a consequence of the parametric construction in Section 4.2, 
triangle areas of M° shrink by a factor of 4 under subdivision, and hence 
the Jacobian of <r t is X A. Each of the entries (l r )hi in * J therefore has one or 
more corresponding entries (I J ' + 1 ) h 'i' i* 1 V\ U P to a factor of V4; that is, Va 
(P') hi = {P f + 1 ) h . r . 
The resulting m X m matrix equation 

p" = (p/)q/+ip/ (11) 

represents a homogeneous system of m 2 equations in the m 2 unknown entries 
of I J '. Due to the symmetry of V \ the system reduces to m(m + l)/2 
homogeneous equations in as many unknowns. A square inhomogeneous 
system is produced once an absolute scale is chosen for the homogeneous 
system. We typically set the scale by requiring that the sum of the entries of I J 
is 1 — this is equivalent to selecting an area form that assigns unit area to M°. 
Although we have been unable to prove that the system is uniquely solvable 
(equivalently, that the system matrix is of full rank), this has been true in 
hundreds of cases we have tried, including all those used to generate the 
figures in this article. We therefore conjecture that under mild conditions on 
the subdivision scheme the system is uniquely solvable. 
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Fig. 6. Computing the inner product (<j> J i(x) f <^ +1 (x)> for piecewise linear scaling functions, 
(a) the graph of <j>i{x) is drawn using solid lines, and the graph of <f>i+i(x) is drawn using 
dashed lines; (b) <M(*> and <M-t-x(*) after refinement. The original problem reduces to that of 
computing the nine inner products between a solid function and a dashed function. 

Once the entries of I J ' have been determined, the remaining inner 
product matrices F'" 1 , F'~ 2 , . . . , 1° can be successively determined via 
Equation 11. 

As a simple example of the integration process, consider the one-dimen- 
sional case where the scaling functions 4^{x) are piecewise linear "hat" 
functions parametrized over the infinite real line, as illustrated in Figure 6. 
For this case, a scaling function with knots on the even integers can be 
refined according to: 



where scaling functions at level j + 1 are hat functions with knots on the 
integers. 

For any level j, these functions lead to only two nonzero cases of inner 
products: 

(1) (<f>i(x) t 4>i(x)) — the inner product of a scaling function with itself. 

(2) (<f>{(x), <l>i+i(x)) — the inner product of a scaling function with its 
immediate neighbor. 

All other possible inner products are either equivalent to one of these, or 
zero, due to the local support of <f>{(x). 
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Choosing our scale such that case 1 evaluates to 1, explicit integration of 
the piecewise linear functions shows that the inner product for case 2 is Va. 
We now determine these same integrals using the general integration 
technique. 

We define the unknown inner product values for case 1 and case 2 to be 

*i : = <M*)> <H(*)> 

(13) 

x* : = <M + i(*)>. 

Figure 6(a) illustrates case 2 — the inner product of two neighboring piece- 
wise linear scaling functions at level j with knots on the even integers. 
The next step is to refine with Equation 12: 

Bilinearity of inner products allows this to be expanded to: 
4 Z 

+ 7 M'-M, 4>%h(x)) 

4 

+ £ i+ZKx), 4%h(x)) + (t&ix), ^h(x)) + \ <*§i( x ), 4>& a (x)) 
2 ^ 

+ 7 to&V*), ^a+V*)) + £ <*'a + 1 i(x), 4>$ti*)> 
4 ^ 

+ 7 (<k J 2iUx), <t> J ^ 3 (x)). 
4 

The result after refinement is shown in Figure 6(b). Through refinability, 
Equation 13 has been expressed in terms of the nine possible inner 
products involving a solid function and a dashed-line function. Most of 
these terms drop out because the supports of their respective functions do 
not overlap: 

x 2 = \ <ti +1 (x), 4>tiKx)) + i (♦/tfu), h:im> + \ (titiM, n:i(x)>. 

(14) 
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The inner products in the first and last terms of Equation 14 are 
equivalent to case 2. However, they are parametrized over a narrower 
domain, and are therefore reduced by a scale of Vfe (recall that the 
appropriate scale for surfaces is V4). Likewise, the inner product of the 
middle term is similar to case 1. We can therefore rewrite Equation 14 in 
terms of the unknowns x x and x 2 : 

1 1 

Similar analysis of case 1 yields: 

3 

*i = 7Xi + x 2 . 
4 

After subtracting out the unknowns on the left, these equations can be 
written as a homogeneous linear system of the form 



8 2 
1 

— 1 
\ 4 / 

Arbitrarily setting x 1 := 1 again yields the inner product for case 2: 
x 2 = Vi. This agrees with the result obtained earlier through explicit 
integration. 

6. MULTI RESOLUTION ANALYSIS BASED ON SUBDIVISION 

In previous sections we have established nested linear spaces and an inner 
product relative to a subdivision rule. We are now in a position to construct 
wavelets, that is, a set of functions ^'(x) = dK(x), <H(x), • • •) that span 
wavelet spaces W J (M°). 

It is of significant practical importance that the analysis and synthesis 
filters associated with these wavelets are constructed and applied in linear 
time. This practical concern drives much of the development in this section, 
the remainder of which is structured as follows. In Section 6.1 we first 
describe a simple construction leading to semiorthogonal wavelets. Al- 
though these wavelets do not satisfy the linear time requirements, a 
variant of the construction is used in Section 6.2 to produce biorthogonal 
wavelets. In Section 6.3 it is shown that for interpolating subdivision 
schemes the biorthogonal wavelets possess linear time analysis and synthe- 
sis procedures. 

6.1 The Semiorthogonal Construction 

Our construction of semiorthogonal wavelets consists of two steps. First, we 
build a basis for V J+l (M°) using the scaling functions * 7 (x) together with 
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the new scaling functions Jf j+1 (x) in V" +1 (Af°). It is straightforward to 
show that *'(x) and tf J+1 (x) together span V J+1 (M°) if and only if the 
matrix O j (encoding the subdivision rule around the old vertices) is 
invertible. Most primal subdivision methods, such as polyhedral subdivi- 
sion and the butterfly scheme, have this property. 1 

Given a function S J+1 (x) in V J+1 (M°) expressed as an expansion in the 
basis (4> y (x) J^ +1 (x)), an approximation in V J (M°) can be obtained by 
restriction to * J (x), that is, by setting to zero the coefficients corresponding 
to N J+1 (x). However, this generally does not produce the best least-squares 
approximation. To ensure the best least-squares approximation after re- 
striction to 4> 7 (x), we may orthogonalize the new basis functions N J+l (x) by 
subtracting out their least-squares projection into V J (M°). Expressed in 
matrix form: 



The coefficients cr 7 are the solution to the linear system formed by taking 
the inner product of each side of Equation 15 with $> j (x\ and using the fact 
that <* J (x), *>(x)> = 0: 



where <F, G) stands for the matrix whose i 9 i'th entry is ((F) i9 (G) r ). The 
matrix (® j (x\ ^(x)) is therefore simply V\ and the matrix <4>-> +1 (x), 
Jt J+1 (x)) is a submatrix of I J+1 consisting of those columns corresponding 
to members of Jf" +1 (x). The system of equations may be solved for the 
coefficients a*. As an example, Figure 2 shows the matrix a 0 for the 
tetrahedron. 

Although this construction produces semiorthogonal wavelets — that is, 
the wavelets M^Cx) are orthogonal to the scaling functions in V J (M°) — they 
are not practically useful because analysis and synthesis both require 
quadratic time. To see why, note that the inner product matrix I J in 
Equation 16 must be inverted. Although 1 J is sparse, its inverse is dense, 
implying that the wavelets are globally supported on M°. The synthesis 
matrices Q- 7 mentioned in Section 2 and more fully described in Section 6.3 
are therefore dense, leading to quadratic time synthesis. The methods used 
in Section 6.3 can be used to show that the analysis matrices A J are also 
dense, making analysis a quadratic time process as well. 

As we show in the next two sections, linear time analysis and synthesis 
can be achieved, at least for interpolating subdivision schemes, by modify- 
ing the above construction to produce biorthogonal wavelets. 



1 One notable exception is Catmull-Clark subdivision for vertices of valence three. However, 
the subdivision rule for such vertices can be easily modified to produce an invertible matrix. 



¥'(x) = Jf' +1 (x)-*'"(x)« / . 



(15) 



<*'"(x), *>(x)>a'"= <<Hx), ^ ;+1 (x)), 

= (P') r <*'" +1 (x), JT' +1 (x)> 



(16) 
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6.2 The Biorthogonal Construction 

linear time synthesis can be achieved by using locally supported wavelets, 
the construction of which is a common problem in the wavelet literature, 
and has been handled in various ways. For instance, Daubechies [1988] 
uses Fourier techniques to develop locally supported fully orthogonal 
wavelets, and Chui [1992] relies on the rich theory of B-splines to develop 
minimally supported semiorthogonal spline wavelets. However, neither of 
these approaches is appropriate for our use: Fourier techniques rely on 
translation and dilation, and B-splines rely on planar topologies. 

We currently do not know of an orthogonal or semiorthogonal construc- 
tion possessing linear time analysis and synthesis for arbitrary topological 
domains. Our approach is to relax semiorthogonality and construct instead 
biorthogonal wavelets. That is, we no longer require the wavelets W j (x) to 
be orthogonal to V J (M°), but to preserve good approximation properties, we 
require the wavelets to be "as orthogonal as possible" subject to the linear 
time requirement. (We should note that a similar approach has been taken 
in independent work by Dahmen et al. [1993].) 

Our biorthogonal construction proceeds by selecting the supports of the 
wavelets a priori. The idea is to "partially orthogonalize" each new scaling 
function <^ +1 (x) in Jf / ' +1 (x) by subtracting off a locally supported least- 
squares best projection of it into V J (M°). Stated another way, we allow 
each column of the matrix a J in Equation 15 to have only a constant 
number of nonzero entries. We then determine the nonzero entries of ot 3 by 
minimizing the L 2 norm of the projection of <t>^(x) into V J {M°). 

More concretely, let <£-J +1 (x) denote a new scaling function, and let $ J t 
denote any subset of the indices of the scaling functions in V J (M°); 
typically ${ consists of the indices of scaling functions in V J (M°) supported 
in some neighborhood of <£-{ +1 (x). We then take the wavelet iH(x) to be 

= *{ +1 (x) + S d7) 

where the coefficients a J r i are such that the L 2 norm of the projection of 
i/rj(x) into V J (M°) is minimized. These coefficients are therefore deter- 
mined by solving the following local linear system: 

2 a{j(4>{-(x)> <«W> = <# +1 (x), <f>i(x)>, (18) 

for all k in 

One way to control the support of the wavelets in a symmetric fashion is 
to select the index set ${ as follows. The new vertex v{ +1 of M J+1 is 
associated with an edge u — v of M j . For instance, the vertex numbered 7 
of M 1 in Figure 2 is associated with the edge 2 - 3 of M °. We take i\ to be 
the vertex indices of M J in the k -discs of u and v. (The k -disc around a 
vertex i; of a triangulation is defined to be the set of all vertices reachable 
from v by following k or fewer edges of the triangulation.) Except for Tables 
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Fig. 7. A polyhedral wavelet centered on a vertex of valence 6. 



II and III, all examples in this article were computed using 2-discs. Figure 
7 is a plot of a polyhedral wavelet computed using 2-discs. 

As the size of the discs increases, the wavelet supports grow, and the 
wavelets have empirically been observed to rapidly approach semi-orthogo- 
nality. 

6.3 A Fitter Bank Algorithm 

< Multiresolution analysis on the infinite real line is based on an assumption 
of spatial invariance — intuitively, every place looks like every other place. 
This means that the analysis and synthesis filters can be represented by a 
convolution kernel, that is, by a sequence of real numbers. This is not the 
case for multiresolution analysis on arbitrary topological domains. The 
filter coefficients in general must vary over the mesh, so the filters are 
represented by (hopefully sparse) matrices. 

The analysis and synthesis filters can be conveniently expressed using 
block matrix equations. For any multiresolution analysis the synthesis 
filters are defined by the relation 



The analysis filters can be used to decompose a surface S J+1 (x) in 



(4>;( X )^( X )) = 3>> +1 (x)(P' QO, 



(19) 



and the analysis filters are obtained from the inverse relation 




(20) 
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V J+1 (M°) given by 

S> +1 (x) = £ uj +1 <M +1 (x) (21) 

i 

into a lower resolution part in V J (M°) plus a detail part in W J (M °) 

s> +1 (x) = S + 2 w<*i(x) 

i 

i 

as follows. Let V be as in Equation 5, and let W- 7 denote the corresponding 
matrix of wavelet coefficients w\. We can rewrite Equation 21 in matrix 
form and substitute the definition of the analysis filters. Thus: 

S' +1 (x) = 3>; +1 (x)V' +1 

= *>(x)A>V> +1 + V J (x)B J V J+l 

therefore, 

V' = A'"V +1 (22) 

W' = B'V' +1 . (23) 

Of course, the analysis filters A 7 " 1 and B^ 1 can now be applied to V to 
yield V-*" 1 , W-*" 1 , etc. A similar argument shows that V +1 can be recov- 
ered from V J and Vf J using the synthesis filters: 

v> +l = p; yj + Q'W. (24) 

We shall now develop more explicit expressions for the analysis and 
synthesis filters. It is again convenient to write 4> J+1 (x) in block form as 

(0> +1 (x)J^" +1 (x)). 

It then follows from Equation 15 that the synthesis filters can be written in 
block form as 

<«■<■ *>-(£ 

where 1 denotes the identity matrix. 

The analysis filters are obtained from Equation 20. (Examples of analysis 
and synthesis matrices for the tetrahedron are shown in Figure 2.) To 
achieve linear time analysis and synthesis, the analysis and synthesis 
matrices must be sparse, having only a constant number of nonzero entries 
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Table I. Time to construct analysis niters for various subdivision methods. 



Method 


Continuity 


Asymptotic run time 


Polyhedral 


c° 


0(n) 


Butterfly 


G l 


0(n) 


Loop 


G l 


Speed of sparse linear equation solution 


Catmull & Clark 


G 1 


Speed of sparse linear equation solution 



in each row. If P y and a J are sparse, then from Equation 25 Q J is also 
sparse. Unfortunately, the analysis filters derived from Equation 20 need 
not be sparse. For interpolating subdivision schemes such as polyhedral 
subdivision and the butterfly scheme, the situation is much improved. Such 
interpolating schemes have the property that O j is the identity matrix. In 
this case the resulting filters are 



A A _ / 1 - a*N j a J \ 
B j ) ~ \ -N j 1 /' 



If P J and a J are sparse, then all four filters are also sparse, leading to 
linear time analysis and synthesis. The situation is less desirable for 
methods related to B-splines, such as Loop's scheme and Catmull-Clark 
surfaces. For these subdivision schemes, the synthesis filters are sparse, 
but the analysis filters are dense. This implies that synthesis is still 
possible in 0(n) time, but that the speed of analysis depends on the time to 
invert the sparse analysis matrix of Equation 20, or to solve the related 
sparse linear system. Whether these methods can be made efficient for 
multiresolution analysis is a topic for future investigation. Table I shows 
what we currently know about the time required to develop analysis filters 
for various subdivision methods. 

The filter bank process specialized to the case of piecewise linear func- 
tions is presented in Section 7.1. 

6.4 Stability 

One important question concerns the stability of the synthesis and analysis 
filters. One measure of the stability of a transformation is its L°° norm. The 
L°° norm for any matrix T, denoted by HTHco, is the maximum sum of the 
absolute values of the elements in any row of T, and indicates the 
maximum scaling effect of T on any vector that it transforms. 

For fully orthogonal wavelet bases over regular grids, the L°° norm is 
bounded by a constant, independent of the number of filtering steps. In this 
section, we briefly analyze the L°° norms of the analysis and synthesis 
filters our locally supported biorthogonal wavelets. 

We first consider repeated application of the synthesis filters. This 
transformation maps the coarse-level scaling coefficients V° and the wave- 
let coefficients W°, . . . , W n_1 into the fine-level scaling coefficients V 71 . 
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Expanding equation 24 repeatedly yields 



/ \ 




1 \ 


n-l 


n-l 


n-l 


n p ; 


v°+ S 


n p' 




r-o 


\j-r J 



The L°° norm of this transformation depends on the L°° norm of the product 
Uj~oP j . We claim that this norm is bounded independent of n. By defini- 
tion, convergence of the subdivision scheme defined by the P j implies that 
the entries of 11* "q 1 P j are also bounded independent of n. Since the 
number of entries per row of n"^ 1 P J is independent of n, the L 00 norm of 
Uj=o P- 7 is bounded independent of n. To conclude, we note that V" is the 
sum of 7i H- 1 transformations with bounded L°° norms. Therefore, the norm 
for the entire transformation is at most O(n). As a common special case, 
when the subdivision rule encoded by P J is a convex combination (such as 
with piecewise linear subdivision or with the Catmull-Clark scheme), the 
L°° norm of any P J is exactly 1, implying that the L 00 norm of their product 
is also 1. 

We next consider repeated application of the analysis filter. This trans- 
formation maps the fine scaling coefficients V" into the coarse scaling 
coefficients V° and the wavelet coefficients W°, . . . , W 1 " 1 . Expanding 
equations 22 and 23 yields 

n-l 

v° = n a> v. 

7 = 0 



rt-1 

W' = B' n A'"'V\ 

The stability of these transformations depends on the L°° norm of n^Jo 1 A J . 
Note that the product of these matrices is the projection operator that maps 
V" into V°. We conjecture that the norm for this transformation is also 
bounded independent of n. 

In practice, the norms for these filters appear to be very small, implying 
good stability. Table II gives the V° norms for the analysis filter A 7 over a 
varying number of levels using a varying size support. These A J are built 
using the wavelet approximations described above, for the case of piecewise 
linear subdivision over the octahedron. In the table, size gives the disc size 
chosen for the support of the wavelet approximations. (The &-disc around a 
vertex v of a triangulation is defined to be the set of all triangles whose 
vertices are reachable from v by following k or fewer edges of the triangu- 
lation.) All examples in this article were generated using wavelet approxi- 
mations supported on 2-discs. 
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Table II. IT Norms for the Analysis Filters A J that Arise from Piecewise Linear 
Subdivision on the Octahedron. 





1-disc 


2-disc 


3-disc 


l|A U ||oo 


1.25 


1.3125 


1.3125 




1.27174 


1.40771 


1.4434 


l|A*||co 


1.27174 


1.40859 


1.45371 


l|A 3 ||oo 


1.27174 


1.40859 


1.45372 


l|A 4 ||co 


1.27174 


1.40859 


1.45372 



Table III gives the cumulative effect of the analysis filters. It gives the 
L°° norms for each stage of the composition of the A- 7 , beginning with A 4 , 
down to A 0 . 

The L°° norm for the synthesis filters P J arising from piecewise linear 
subdivision is simply 1 for every level, and is therefore not shown in a 
table. 

Note that these results on stability are only numerical. Necessary and 
sufficient conditions guaranteeing stability in the irregular case are very 
difficult to establish. Schroder and Sweldens [1995] suggest constraining 
the resulting wavelets to have a zero order moment (i.e., the integral of the 
wavelet is zero). In the regular univariate case, this condition is known to 
be necessary for stability [Daubechies 1992]. However, it is unknown 
whether this is a necessary condition for stability in the irregular case. 

If a zero order moment is necessary for stability, enforcing such a 
constraint is easy: During the partial orthogonalization process, one simply 
adds a linear constraint on the unknown coefficients a{ r of the wavelet 
that forces the wavelet to have a vanishing integral. For more information, 
Dahmen [1994] gives a much deeper treatment of the stability of multireso- 
lution schemes over irregular grids, and derives some fundamental approx- 
imation-theoretic results. 

7. WAVELET COMPRESSION OF SURFACES 

Multiresolution analysis is widely used for data compression applications. 
Mallat [1989] and DeVore et al. [1992], among others, use wavelet tech- 
niques for efficient image compression. Finkelstein and Salesin [1994] 
develop a wavelet-based method for approximating B-spline curves within 
an L 00 tolerance. Other examples of wavelet-based compression techniques 
abound [Berman et al. 1994; Chui and Shi 1992; Lucier 1992; Mallat and 
Hwang 1992]. 

Using wavelet techniques, lossy compression of a function is typically 
implemented with a three-stage algorithm: 

(1) Filter bank decomposition. The original function, represented by the 
sequence v J > is decomposed into a coarse-level approximation v° to- 
gether with wavelet coefficients w° 9 . . . , w j ~ x at the levels from 0 to 
j - 1. 
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Table HI. IT Norms for Successive Products of the A J \ Beginning at A 4 



Matrix product 


1-disc 


2-disc 


3-disc 


A 4 


1.27174 


1.40859 


1.45372 


A 3 * A 4 


1.61574 


1.95076 


2.04751 


A** A 3 * A 4 


2.02336 


2.50179 


2.60807 


A 1 * A 2 * A 3 * A 4 


2.35466 


2.78397 


2.86811 


A°*A 1 *A 2 *A 3 *A 4 


2.22687 


2.63381 


2.75473 



(2) Selection, A subset d of the detail is chosen from each sequence 
w°, . . . , u;- 7-1 , according to some measure of its importance. 

(3) Reconstruction, The selected detail d, in the form of scaled wavelets at 
various levels, is added back to the coarse-level approximation. 

Throughout the rest of this section, we will discuss the application of 
subdivision wavelets to the compression of complex surface data. The 
common special case of polyhedral decomposition is treated in Section 7.1. 
Next, Section 7.2 discusses various rules for selecting wavelet coefficients. 
Section 7.4 discusses the preprocessing step sometimes necessary for 
converting general input into a form appropriate for compression. Finally, 
Section 7.5 illustrates the results of compression for two complex polyhe- 
dral data sets. 

In this section, we exclusively consider the approximation of functions 
defined over triangular meshes. Quadrilateral methods require separate 
algorithms which are nevertheless similar in flavor. 

7.1 Polyhedral Implementation 

Many applications of subdivision wavelets involve the special case of 
piecewise linear subdivision — an example is polyhedral compression. The 
techniques outlined in Section 4 apply to polyhedra, but it is possible to 
exploit the simplicity of piecewise linear subdivision to derive an even more 
efficient implementation. Moreover, decomposition and reconstruction may 
be achieved in linear time for polyhedral surfaces without the need for a 
sparse matrix representation. (However, it is still necessary to solve a 
small linear system for the local neighborhood of wavelet coefficients.) 

The chief property of piecewise linear subdivision that leads to simpler 
algorithms is tight locality. Unlike more complex subdivision rules, the 
support of a hat function built around a vertex v does not extend beyond v*s 
neighboring vertices. In this section, we show that this leads to a simpler 
technique for generating the filters A and B. 

As is discussed in Section 6.3, the wavelet coefficients Vf J may be found 
through the matrix product H-ty 7 '*** 1 . Construction of B- 7 is greatly simpli- 
fied for polyhedral subdivision, where the wavelet coefficient w associated 
with a vertex v (where v is fine-level vertex added on the edge pq) may be 
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Fig. 8. Determining the wavelet coefficient w around vertex v, with parents p and q. The 
midpoint between p and q is m, and the resulting wavelet coefficient is the difference w 
between m and v. 

derived by the simple rule: 



This equation is shown geometrically in Figure 8. It leads to a very simple 
algorithm for computing the wavelet coefficients — without the need for 
implementing sparse matrix multiplication. 

Determining all level j - 1 wavelet coefficients w J ~ x using Equation 26 
is the first half of the level j to level j - 1 decomposition process. It is still 
necessary to derive the filter A* 7-1 that gives the scaling functions for the 
approximation at level j - 1. This is not as easy, but one may still take 
advantage of the simplicity of the piecewise linear representation in order 
to achieve a procedure that is more efficient than the general case. 

Once the wavelet coefficients w-^' 1 are derived, the coarser-level approx- 
imation at level j — 1 may be determined by subtracting from the level j 
function the effect of every scaled wavelet term i^ -1 . The problem 

then reduces to determining the wavelets iH" 1 * 

The ma y be computed using the techniques of Section 4, but 

specifically adapted for linear subdivision. In particular, the inner product 
mask for the linear case is shown in Figure 4. This mask is especially easy 
to determine for polyhedra, because the central value is simply the valence 
of the vertex. 

7.2 Coefficient Selection 

Once the surface has been decomposed via the filter bank, the next step in 
surface compression is to select appropriate coefficients to add back to the 
base mesh. 

Before examining the wavelet coefficients for selection, it is useful to first 
normalize them. Normalization allows coefficients at differing subdivision 



1 

w = v - - (p + <?). 



(26) 
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levels to be equitably compared, despite the different domain sizes over 
which they are defined. 

Least-squares (L 2 ) normalization of the wavelet coefficients assigns a 
weight to a wavelet coefficient that indicates how much least-squares error 
occurs when the coefficient is left out of the reconstruction. A wavelet tH can 
be normalized to a "unit length" wavelet i# with the following equation: 

<M 

When the wavelet ty\ has an associated wavelet coefficient w J h the L 2 
normalized coefficient w{ is therefore assigned the value 

All examples presented in this article use selection based upon this L 2 
normalization of the wavelet coefficients. 

There are several possible techniques for selecting coefficients. These 
include: 

— Threshold testing. Perhaps the simplest selection technique is to simply 
choose all wavelet coefficients whose magnitude is greater than some 
coefficient e [Donoho 1994]. Although thresholding is quite simple, the 
results are of remarkable quality, as we will see in Section 7.5. In the 
reconstruction, areas of high curvature are sampled more densely than 
relatively flat regions. 

— L 2 progressive refinement. In certain applications, it is useful to ensure 
that the most important information is reconstructed first. When the L 2 
normalization described above is used, it can be shown that simply 
supplying the c largest coefficients in decreasing order is sufficient to 
produce a sequence of approximations, each of which is the best possible 
least-squares approximation using only c coefficients [Donoho 1994]. The 
wavelet coefficients must first be sorted by their magnitude, which 
implies an 0(n log n) run time in the input size. 

— Maximum error: L°° reconstruction. An approximation constructed ac- 
cording to the L 2 norm may still contain an arbitrarily large error in a 
sufficiently small region. An alternative approach uses the L°° norm to 
guarantee that no part of the reconstruction is farther than a user- 
defined tolerance from its corresponding point on the original input. 

L°° reconstruction begins with the complete, unreduced reconstruction, 
and attempts to remove successive wavelet coefficients until no more are 
possible. L°° wavelet compression of complex polyhedral surfaces is 
further studied in Eck et al. [1995]. 

—Location. In some applications, it may be useful to compress only a 
portion of a model while the rest of it is outside the user's view. For 
example, a user may be interested in viewing a specific location on a map 
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of the Earth, but uninterested in viewing geography elsewhere. In this 
case, the location information associated with a wavelet makes it possible 
to apply wavelet compression to the region of interest, but to avoid 
processing irrelevant regions of the model. 

Selection by location is usually applied in conjunction with another 
selection method. Only those coefficients in the relevant region are 
candidates for the secondary selection technique. 

7.3 Reconstruction Algorithms 

To review, the result of decomposition is the base mesh and the coefficients 
of the wavelets at various levels of subdivision. The information at the 
vertices of the base mesh stores values of a very few scaling functions, 
providing a very coarse least-squares approximation of the input function. 
The wavelets are functions that represent the missing detail at each 
finer-level vertex. In order to build a wavelet approximation of the original 
input, we can add back to the base mesh a selected subset of these 
finer-level functions. 

Reconstruction of the approximation is done by adding back to the base 
surface the scaled wavelet associated with each selected coefficient. Details 
of the reconstruction phase are given explicitly in Lounsbery [1994]. 

The effect of reconstruction is to produce vertices on the approximation 
surface that represent the addition of the selected wavelets to the base 
surface. However, this form is not appropriate for many applications that 
require an actual surface representation. If the approximation is to be 
rendered as a polyhedral surface, it is necessary to generate polygons over 
the approximated surface. A separate top-down recursive triangulation 
algorithm is required in order to connect the points on the reconstruction 
into polygonal form. 

When c is the number of selected coefficients and j is the subdivision 
depth of the input function, the asymptotic time for reconstruction and 
triangulation is at worst 0(j c), but in practice usually much closer to 
0(c). These values are borne out by the empirical evidence of run times for 
reconstruction. 

7.4 Subdivision Connectivity 

The compression techniques detailed in this section are useful for applica- 
tions requiring decomposition of a function in V J (M°), where M° may be a 
mesh of any topological type. An implicit assumption is that the connectiv- 
ity of the input mesh must have the form of a mesh M j that results from 
subdividing a simple mesh M° j times. We call this property subdivision 
connectivity. 

For many applications, such as global illumination or surface editing, 
producing input with subdivision connectivity is a simpler matter. For 
other purposes, including multiresolution compression of an arbitrary 
mesh, an initial preprocessing step is required to convert the input into an 
approximation with the necessary connectivity. Such a conversion algo- 
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rithm is presented in Eck et al. [1995] for the general case. The examples 
presented in this article were converted with an algorithm written for the 
special case of converting input taken in cylindrical sections into points on 
a subdivided octahedron. 

7.5 Polyhedral Examples 

In this section, subdivision wavelets are applied to two compression prob- 
lems: the compression of a polyhedral model consisting of over 32,000 
triangles, and compression of a piecewise linear representation of a color 
function defined on over one million points on the globe. The data for these 
examples were resampled from cylindrical sections onto a subdivided 
octahedron. 

7.5.1 Geometric Data. The input for the first example (shown in Figure 
9(a)) is a polyhedral mesh consisting of 32,768 triangles whose vertices 
were resampled from laser range data originally provided through the 
courtesy of Cyberware, Inc. The triangulation was created by recursively 
subdividing an octahedron six times. The octahedron therefore serves as 
the domain mesh M°, with the input triangulation considered as a para- 
metric function S(x), x G M° lying in V^Af 0 ). More precisely, if v 6 t denotes 
the vertices of the input mesh, S(x) can be written as 

S(x) = 3> 6 (x)V 6 xGM° 

where the scaling functions 4> 6 (x) are the (piecewise linear) functions 
defined through polyhedral subdivision. 

The locally supported wavelet approximations ^{(x) for this example are 
chosen to be supported on 2-discs. The filter bank process outlined in 
Section 6.3 can be applied in linear time to rewrite S(x) in the form 

5 

5(x) = a>°(x)v° + £ *'"(x)W. 

The first term describes a base shape as the projection of S(x) into V°(M°), 
which in this case is an approximating octahedron with vertex positions 
given by the six rows of V°. For this data, the decomposition stage of the 
filter bank runs in about 14 seconds on an SGI Indigo 2 Extreme. 

Approximations to the original mesh S(x) can be easily obtained from the 
wavelet expansion using coefficients selected according to the threshold 
rule (see Section 7.2). The models in Figure 10(b), (d), and (f) are com- 
pressed to 1%, 13%, and 32%, respectively. Notice that thresholding causes 
the mesh to refine in areas of high detail, while leaving large triangles in 
areas of relatively low detail. An implementation of the entire decomposi- 
tion and reconstruction process that produces Figure 10(d) runs in about 53 
seconds from input to output. 
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(a) Distant view of mesh in (b). (b) 105 wavelets, 240 triangles. 




(c) Mid-range view of mesh in (d). (d) 1,966 wavelets, 4,272 triangles. 




(e) Close-up view of mesh in (0. (0 5,054 wavelets, 10.704 triangles. 

Fig. 10. Wavelet approximations of the Spock polyhedron. Left: Gouraud-shaded wavelet 
approximations. Right: Flat-shaded close-ups showing structure of approximations. 
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(a) Full-resolution model. (b) Mesh for (a). 



Pig. 11. Close-up of Earth data before compression. 

Figure 10 also illustrates the use of wavelet approximations for auto- 
matic level-of-detail control in rendering. The original Spock polyhedron is 
shown in full resolution in Figure 9(a). Views of the full-resolution model 
from various distances are shown in the rest of Figure 9. When viewing the 
input polyhedron at these distances, it is inefficient and unnecessary to 
render all 32,000 triangles. The approximations shown in the left column of 
Figure 10 may instead be used without significantly degrading image 
quality. 

Suddenly switching between models of different detail in an animation 
often produces a noticeable jump. This problem is easily mended by using a 
wavelet expansion where the wavelet coefficients are treated as continuous 
functions of the viewing distance. This simple technique allows the object 
geometry to smoothly change its appearance as the viewing distance 
changes. This method has proved successful shown in a frame-by-frame 
animation we have produced. 

7.5.2 Color Data. Subdivision wavelets may be applied to more general 
functions over surfaces than geometric data. Figure 12 demonstrates 
another application — that of compressing a function on the sphere. In this 
example, elevation and bathymetry data obtained from the U.S. National 
Geophysical Data Center was used to create a piecewise linear pseudocol- 
oring of the sphere. The resulting color function was represented by 
2,097,152 triangles and 1,048,578 vertices. The full-resolution pseudocolor- 
ing was too large to be rendered on an SGI Indigo 2 Extreme with 128 MB, 
and is therefore not shown in its entirety in Figure 12. Instead, Figure 11 
shows a close-up of a region that is compressed in Figure 12. An apprecia- 
tion for the density of the data can be obtained from Figure 11(b), where 
even at close range the mesh lines of the original uncompressed data are so 
dense that the image appears almost completely black. 

The approximations shown in Figure 12(a)-(f) were produced using 
subdivision wavelet compression. Figure 12(a) shows a distant view of the 
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(a) Distant view of (b). (b) Compression to 0. 1%. 




(e) Close-up view of (0- (0 Compression to 16%. 

Fig. 12. Approximating color as a function over the sphere. 
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(a) The original smooth surface. (b) Its wavelet compression to 1 6%. 




(c) A coarse-level change of (a). (d) A finer-level modification. 



Fig. 13. Compression and multiresolution editing of a smooth surface. 

Earth using an approximation of only 0.1% (the mesh is shown in Figure 
12(b)). Likewise, Figures 12(c) and (d) show the result of compression to 2% 
for a medium-range view. At close range the compression to 16% in (e) is 
nearly indistinguishable from the full-resolution model in Figure 11(a). A 
comparison of the meshes shown in Figure 11(b) and Figure 12(f) reveals 
the striking degree of compression achieved in this case. 



8. SMOOTH SURFACE MODELING 

This section describes how subdivision wavelets can be used to compress 
and edit smooth subdivision surfaces. Although our examples concentrate 
on the modeling of smooth surfaces based on the Dyn et al. [1990] butterfly 
scheme, the techniques described in this section may be applied to general 
subdivision surfaces, at the cost of potentially quadratic processing time in 
the case of Loop's method. 
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In Section 7, we saw examples of polyhedral compression using wavelets. 
In this section, we examine compression of surfaces defined by smooth 
subdivision rules. 

8.1 Smooth Surface Compression 

The result of smooth surface compression is shown in Figure 13. The 
surface in Figure 13(a) shows the level 6 Spock data set from Section 7.5 
after it has been subdivided using the tangent-plane smooth butterfly 
subdivision rule. Figure 13(b) shows the result of compression using 
wavelets that were developed from the butterfly scheme. The surface of 
Figure 13(a) has been compressed to a representation (depicted in Figure 
13(b)) that may be stored using only 16% of the original coefficients. 

8.2 Multiresolution Editing 

Subdivision wavelets may also be used for editing shapes at multiple 
resolutions. Although we have not fully implemented multiresolution edit- 
ing, Figures 13(cMd) show an example of editing the smooth shape seen in 
Figure 13(a). 

Figure 13(c) shows the effect of changing a single scaling function 
coefficient on the level 0 base octahedron. Because finer-level vertices in 
the same region are defined relative to the coarser shape, they move along 
with the modification. However, the geometry in areas away from the front 
of the bust is not affected. 

It is also possible to locally modify the shape at a finer level by changing 
the value of a wavelet coefficient at that level. The result of modifying a 
single level 3 wavelet coefficient is shown in Figure 13(d). 

The changes shown in Figure 13(c)-(d) were created by simply modifying 
. ^ a single value in the wavelet representation. It is preferable to make these 
changes under a more powerful user interface. Finkelstein and Salesin 
[1994] discuss interfaces that are more appropriate for wavelet editing of 
curves and surfaces. 

9. SUMMARY AND FUTURE WORK 

In this article, we have established a theoretical basis for applying mul- 
tiresolution analysis to surfaces of arbitrary topological type. These results 
hold for any local, uniformly convergent, continuous, primal subdivision 
scheme, including polyhedral subdivision, the butterfly scheme [Dyn et al. 
1990], Loop's scheme [1987], and Catmull-Clark surfaces [1978]. The 
results also hold for piecewise smooth subdivision as described in Hoppe 
[1994] and Hoppe et al. [1994], and for open surfaces possessing bound- 
aries. 

There are numerous areas for future work: 

— Linear time sparse matrix inversion. As explained in Section 6.3, linear- 
time reconstruction for general subdivision rules depends upon linear- 
time solution of general sparse linear systems. Hence, finding an 0(n) 
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sparse linear system solver is of value for efficiently implementing 
subdivision wavelets for more general subdivision schemes, such as the 
G 1 methods of Loop and Catmull-Clark. 
— Semiorthogonal wavelets. The wavelets developed in Section 6 are 
biorthogonal and linear-time for interpolating subdivision. It would be 
interesting to determine if linear-time semiorthogonal wavelets also 
exist. 

— Simplifying the topological type. The surface decomposition developed 
herein retains the topological type of the input surface. When the input is 
a relatively simple object with many small holes, it is more often 
desirable to decompose the input into a "topologically simpler" surface; 
that is, one with lower genus or fewer boundary curves. 
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